What if we wanted to look at population in 20 years given an initial condition
Two options
explicit solution to differential equation is known; e.g. you can integrate both sides of the equation! Not always possible but lets look at a case where it is possible
must be solved by iteration; this is what we do when we can’t integrate both sides
## function (T, P0, r, K)
## {
## P = P0 * exp(r * T)
## if (P > K) {
## P = K
## }
## return(P)
## }
# gives population after any time given an initial population
# 20 rabbits, growth rate of 0.01 how many in 30 years
exppop(T=30, P0=20, r=0.01, K=1000)## [1] 26.99718
# if we want to see how population evolves over time - generate a time series by running our model for each point in time
initialrabbits = 20
years = seq(from=1, to=100, by=2)
Ptime = years %>% map_dbl(~exppop( P0=initialrabbits, r=0.01, K=1000, T=.x))
# keep track of what times we ran
Ptime = data.frame(P=Ptime, years=years)
ggplot(Ptime, aes(years,P))+geom_point()+labs(x="years",y="Rabbit Population")Using a solver….when you can’t do the integration by hand :)
For example, if you added a carrying capacity as threshold where it stops growing
In this case we integrate by iteration - approximates analytic integration
Implement the differential equation as a function that
*inputs time, the variable(s) and a parameter list
(it needs time even though you don’t use it)
My convention: name derivative functions starting with d to remind myself that they are computing a derivative
Only works for Ordinary Differential Equations - single independent variable (in our case time)
Partial differential equations - more than 1 independent variable (e.g x and y if changing in space)
R has a solver called ODE for solving ordinary differential equations from package desolve
## function (time, P, r)
## {
## dexpop = r * P
## return(list(dexpop))
## }
library(deSolve)
initialrabbits = 20
years = seq(from=1, to=100, by=2)
# run the solver
Ptime = ode(y=initialrabbits, times=years, func=dexppop,parms=c(0.01))
head(Ptime)## time 1
## [1,] 1 20.00000
## [2,] 3 20.40404
## [3,] 5 20.81623
## [4,] 7 21.23675
## [5,] 9 21.66576
## [6,] 11 22.10344
colnames(Ptime)=c("year","P")
# notice that there are additional pieces of information year, including the method used for integration
attributes(Ptime)## $dim
## [1] 50 2
##
## $dimnames
## $dimnames[[1]]
## NULL
##
## $dimnames[[2]]
## [1] "year" "P"
##
##
## $istate
## [1] 2 52 105 NA 6 6 0 36 21 NA NA NA NA 0 1 1 NA NA NA
## [20] NA NA
##
## $rstate
## [1] 2.00000 2.00000 99.98839 0.00000 1.00000
##
## $lengthvar
## [1] 1
##
## $class
## [1] "deSolve" "matrix"
##
## $type
## [1] "lsoda"
# this also means you need to extract just the data frame for plotting
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")# this also works (of course function can be by order)
Ptime=ode(initialrabbits, years, dexppop,0.03)
colnames(Ptime)=c("year","P")
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")You can play a bit with changing your function to something that you can’t integrate “by hand”
BUT we might want more parameters
to work with ODE, parameters must all be input as a single list; similar to how we return multiple outputs from a function
see example below..lets add a carrying capacity
## function (time, P, r)
## {
## dexpop = r * P
## return(list(dexpop))
## }
## <bytecode: 0x7fa14a0ce3c8>
# create parameter list
initalrabbits=2
newparms = list(r=0.03, carry_capacity=300)
#apply solver
results=ode(initialrabbits, years, dexppop_play,newparms)
head(results)## time 1
## [1,] 1 20.00000
## [2,] 3 21.23675
## [3,] 5 22.54993
## [4,] 7 23.94435
## [5,] 9 25.42499
## [6,] 11 26.99719
# add more meaningful names
colnames(results)=c("year","P")
# plot
ggplot(as.data.frame(results), aes(year,P))+geom_point()+labs(y="Population", "years")# try again with different parameters
alternativeparms = list(r=0.04, carry_capacity=500)
results2=ode(initialrabbits, years, dexppop_play,alternativeparms)## DLSODA- At current T (=R1), MXSTEP (=I1) steps
## taken on this call before reaching TOUT
## In above message, I1 = 5000
##
## In above message, R1 = 81.5108
##
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
## time 1
## [1,] 1 20.00000
## [2,] 3 21.66575
## [3,] 5 23.47022
## [4,] 7 25.42499
## [5,] 9 27.54257
## [6,] 11 29.83650
colnames(results2)=c("year","P_parm2")
# plot
ggplot(as.data.frame(results2), aes(year,P_parm2))+geom_point()+labs(y="Population", "years")# compare by combining into a single data frame
both = inner_join(as.data.frame(results), as.data.frame(results2))## Joining, by = "year"
both_p = both %>% gather(key="model", value="P", -year)
ggplot(both_p, aes(year,P, col=model))+geom_point()+labs(y="Population", "years")What is ODE doing? (iterating in ‘smart ways’)
Similar to “difference equations”
Population models can be discrete (rather than continuous) So we could implement them as difference equations and iterate
source("../R/discrete_logistic_pop.R")
# notice how a for loop is used to iterate
# how many rabbits after 50 years given a growth of 0.1
# starting with 1 rabbit - but a carrying capcity of 500
discrete_logistic_pop## function (P0, r, K, T = 10)
## {
## pop = P0
## for (i in 1:T) {
## pop = pop + r * pop
## pop = ifelse(pop > K, K, pop)
## }
## return(pop)
## }
## [1] 11.4674
Remember we have 3 ways now to calculate population
analytical solution - based on integration (exppop.R) BEST
using an ode solver for numerical approximation (exppop_play.R)
numerical integration using in discrete steps (discrete_logistic_pop.R)
## [1] 12.18249
## [1] 12.18249
## [1] 11.4674
# why are they different
# look at trajectories
growth_result = data.frame(time=seq(from=1,to=100))
growth_result$Panalytic = growth_result$time %>% map_dbl(~exppop( P0=1,r=0.05, K=200,T=.x ))
growth_result$Pdiscrete = growth_result$time %>% map_dbl(~discrete_logistic_pop( P0=1,r=0.05, K=200,T=.x ))
tmp = growth_result %>% gather(key="Ptype",value="P",-time)
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()Finally look at continuous derivative using ODE solve Needs initial condtions differential equation *parameters
## function (time, P, parms)
## {
## dexpop = parms$r * P
## dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
## return(list(dexpop))
## }
# set up using the same parameters
pcompare = list(r=r,carry_capacity=K)
# now run our ODE solver
result = ode(P0, growth_result$time, dexppop_play, pcompare)
head(result)## time 1
## [1,] 1 1.000000
## [2,] 2 1.051273
## [3,] 3 1.105172
## [4,] 4 1.161835
## [5,] 5 1.221404
## [6,] 6 1.284027
# we already have time - so just extract population
growth_result$Pdifferential = result[,2]
# comapre all 3 approaches
tmp = growth_result %>% gather(key="Ptype",value="P",-time)
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()All differential and difference equations are approximations The analytical solution is exact
Notice that differential equations is a bit more accurate!
Diffusion can be implemented as a partial differential equation Complicated to solve - but there are tool in R for specific types of partial differential equations
More info on differential equations in R
But we can appoximate this with a difference equation - and iterative to get an estimate of how diffusion works (see class lecture pdf)
Now for partial derivatives - time and space! it gets much more tricky …beyond this course
Example of Diffusion - difference equation implementation to see what some issues can be
source("../R/diffusion.R")
# run our diffusion model (iterative difference equation) with initial concentration of 10, for 8 timestep (size 1m), and 10 space steps (size 1s)
# using diffusion parameters 0.5 s/m2, 10 m2
result = diff1(initialC=10, nx=10, dx=1, nt=8, dt=1, D=0.5, area=10)
# a list is returned with our 3 data frames for concentration (conc), qin and qout
result## $conc
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 10.000000 0.000000 0.000000 0.0000000 0.0000000 0.000000000 0.000000000
## [2,] 7.500000 2.500000 0.000000 0.0000000 0.0000000 0.000000000 0.000000000
## [3,] 6.250000 3.125000 0.625000 0.0000000 0.0000000 0.000000000 0.000000000
## [4,] 5.468750 3.281250 1.093750 0.1562500 0.0000000 0.000000000 0.000000000
## [5,] 4.921875 3.281250 1.406250 0.3515625 0.0390625 0.000000000 0.000000000
## [6,] 4.511719 3.222656 1.611328 0.5371094 0.1074219 0.009765625 0.000000000
## [7,] 4.189453 3.142090 1.745605 0.6982422 0.1904297 0.031738281 0.002441406
## [8,] 3.927612 3.054810 1.832886 0.8331299 0.2777100 0.064086914 0.009155273
## [,8] [,9] [,10]
## [1,] 0.0000000000 0 0
## [2,] 0.0000000000 0 0
## [3,] 0.0000000000 0 0
## [4,] 0.0000000000 0 0
## [5,] 0.0000000000 0 0
## [6,] 0.0000000000 0 0
## [7,] 0.0000000000 0 0
## [8,] 0.0006103516 0 0
##
## $qout
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 25.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [2,] 12.500000 6.250000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [3,] 7.812500 6.250000 1.562500 0.000000 0.00000000 0.00000000 0.000000000
## [4,] 5.468750 5.468750 2.343750 0.390625 0.00000000 0.00000000 0.000000000
## [5,] 4.101562 4.687500 2.636719 0.781250 0.09765625 0.00000000 0.000000000
## [6,] 3.222656 4.028320 2.685547 1.074219 0.24414062 0.02441406 0.000000000
## [7,] 2.618408 3.491211 2.618408 1.269531 0.39672852 0.07324219 0.006103516
## [8,] 0.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [,8] [,9] [,10]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 0
## [5,] 0 0 0
## [6,] 0 0 0
## [7,] 0 0 0
## [8,] 0 0 0
##
## $qin
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 25.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000
## [2,] 0 12.500000 6.250000 0.000000 0.000000 0.00000000 0.00000000
## [3,] 0 7.812500 6.250000 1.562500 0.000000 0.00000000 0.00000000
## [4,] 0 5.468750 5.468750 2.343750 0.390625 0.00000000 0.00000000
## [5,] 0 4.101562 4.687500 2.636719 0.781250 0.09765625 0.00000000
## [6,] 0 3.222656 4.028320 2.685547 1.074219 0.24414062 0.02441406
## [7,] 0 2.618408 3.491211 2.618408 1.269531 0.39672852 0.07324219
## [8,] 0 0.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000
## [,8] [,9] [,10]
## [1,] 0.000000000 0 0
## [2,] 0.000000000 0 0
## [3,] 0.000000000 0 0
## [4,] 0.000000000 0 0
## [5,] 0.000000000 0 0
## [6,] 0.000000000 0 0
## [7,] 0.006103516 0 0
## [8,] 0.000000000 0 0
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 10.000000 0.000000 0.000000 0.0000000 0.0000000 0.000000000 0 0 0
## [2,] 7.500000 2.500000 0.000000 0.0000000 0.0000000 0.000000000 0 0 0
## [3,] 6.250000 3.125000 0.625000 0.0000000 0.0000000 0.000000000 0 0 0
## [4,] 5.468750 3.281250 1.093750 0.1562500 0.0000000 0.000000000 0 0 0
## [5,] 4.921875 3.281250 1.406250 0.3515625 0.0390625 0.000000000 0 0 0
## [6,] 4.511719 3.222656 1.611328 0.5371094 0.1074219 0.009765625 0 0 0
## [,10]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0
# or if you prefer this orientation (Distance on x axis)
filled.contour(t(result$conc), ylab="Time", xlab="Distance")# changes diffusivity and other parameters particularly
# diffusivity, dx and dt
res=diff1(initialC=10,nx=10,dx=1,nt=10,dt=30,D=0.006,area=1)
filled.contour(res$conc, xlab="Time", ylab="Distance")# we can also see how much material moved from place to place each time step
filled.contour(res$qin, xlab="Time", ylab="Distance")# what if we increase diffusivity
resfast=diff1(initialC=10,nx=10,dx=0.5,nt=10,dt=10,D=0.08,area=1)
filled.contour(resfast$conc, xlab="Time", ylab="Distance")# Discretization Issue Example
resunstable=diff1(initialC=10,nx=10,dx=1,nt=10,dt=10,D=0.8,area=1)
filled.contour(resunstable$conc, xlab="Time (fraction of hour)",ylab="Distance Along Path (m)", main="Pollutant Diffusion")# this illustrates the problem with difference equations (and the challenges that methods for numerical integration try to overcome)
# if things are changing quickly we need to use much smaller time, space steps to avoid overshoot and instability
# so lets cut our step size by 10 (dt) (but then add 10 more steps (nx) to cover the same distance)
resunstable=diff1(initialC=10,nx=100,dx=1,nt=10,dt=1,D=0.8,area=1)
filled.contour(resunstable$conc, xlab="time",ylab="Distance Along Path")